<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta name="generator" content="MATLAB R2016a"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><title>Modelling coupled miscible fluid flow and transport processes </title><style type="text/css">
* {margin: 0; padding: 0;}
body {text-align: start; line-height: 17.2339992523193px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Consolas, Inconsolata, Menlo, monospace; font-style: normal; font-size: 14px; font-weight: normal; text-decoration: none; white-space: normal; }
h1, h2 {font-weight: normal;}
.content { padding: 30px; }

.S0 { margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S1 { line-height: 26.3999996185303px; min-height: 24px; white-space: pre-wrap; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-size: 22px; white-space: pre-wrap; margin-left: 4px; margin-top: 3px; margin-bottom: 15px; margin-right: 10px;  }
.S2 { min-height: 0px; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S3 { line-height: 21px; min-height: 17px; white-space: pre-wrap; font-family: Helvetica, Arial, sans-serif; white-space: pre-wrap; margin-left: 4px; margin-top: 2px; margin-bottom: 9px; margin-right: 10px;  }
.S4 { font-family: Helvetica, Arial, sans-serif; margin-left: 0px; margin-top: 10px; margin-bottom: 20px; margin-right: 0px;  }
.S5 { text-align: left; line-height: 21px; white-space: pre-wrap; white-space: pre-wrap; margin-left: 56px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S6 { line-height: 20.576000213623px; min-height: 20px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-size: 16px; font-weight: bold; white-space: pre-wrap; margin-left: 4px; margin-top: 3px; margin-bottom: 9px; margin-right: 10px;  }
.S7 { margin-left: 3px; margin-top: 10px; margin-bottom: 10px; margin-right: 3px;  }
.S8 { line-height: 15.5926666259766px; min-height: 18px; white-space: nowrap; font-size: 12.6666669845581px; white-space: nowrap; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S9 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; color: rgb(34, 139, 34); font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 45px;  }
.S10 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 45px;  }
.S11 { line-height: 21px; min-height: 17px; white-space: pre-wrap; font-family: Helvetica, Arial, sans-serif; white-space: pre-wrap; margin-left: 4px; margin-top: 10px; margin-bottom: 9px; margin-right: 10px;  }
.S12 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S13 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; color: rgb(160, 32, 240); font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S14 { color: rgb(64, 64, 64); margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S15 { margin-left: 3px; margin-top: 10px; margin-bottom: 4px; margin-right: 3px;  }
.S16 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; color: rgb(0, 0, 255); font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S17 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; color: rgb(0, 0, 255); font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 45px;  }
.S18 { text-align: center; line-height: 21px; min-height: 17px; white-space: pre-wrap; font-family: Helvetica, Arial, sans-serif; white-space: pre-wrap; margin-left: 4px; margin-top: 2px; margin-bottom: 9px; margin-right: 10px;  }

.LineNodeBlock {margin: 10px 0 10px 0;}
.LineNodeBlock+.paragraphNode {margin-top: 10px;}
.lineNode {padding-left: 10px; background-color: #F7F7F7; border-left: 1px solid #E9E9E9; border-right: 1px solid #E9E9E9;}
.inlineWrapper:first-child .lineNode,.inlineWrapper.outputs+.inlineWrapper .lineNode {padding-top: 5px; border-top: 1px solid #E9E9E9;}
.inlineWrapper:last-child .lineNode,.inlineWrapper.outputs .lineNode {padding-bottom: 5px; border-bottom: 1px solid #E9E9E9;}
.lineNode .textBox {white-space: pre;}
.outputGroup {margin: 2px 0 2px 0; padding: 2px 2px 2px 4px;}
.outputRegion {}
.outputParagraph {color: rgba(64, 64, 64, 1); padding: 10px 0 6px 17px; background: white; overflow-x: hidden;}
.inlineWrapper:last-child .outputParagraph {border-bottom-left-radius: 4px; border-bottom-right-radius: 4px;}
.outputParagraph:empty {margin: 0;}
.inlineElement .symbolicElement {}
.embeddedOutputsSymbolicElement .MathEquation {margin-top: 5px; margin-bottom: 5px;}
.embeddedOutputsSymbolicElement .MathEquation.displaySymbolicElement {margin-left: 15px;}
.embeddedOutputsSymbolicElement .MathEquation.inlineSymbolicElement {margin-left: 5px;}
.embeddedOutputsSymbolicElement {overflow-x: auto; overflow-y: hidden;}
.embeddedOutputsSymbolicElement { overflow: initial !important;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px; white-space: pre; word-wrap: initial; min-height: 18px; max-height: 250px; overflow: auto;}
.textElement {padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement { max-height: none !important; overflow: initial !important;}
.veSpecifier {}
.veContainer {}
.veSpecifierBox {height: 400px; width: 500px;}
.veSpecifier .veTable {padding-top: 3px; padding-bottom: 4px;}
.veSpecifierBox .veSpecifier .veContainer {position: relative; width: 100%; height: 370px;}
.veSpecifierBox .dijitDialogPaneContent {width: 97% !important; height: 88% !important;}
.veSpecifier .veTable .rowHeadersWrapper {padding-bottom: 0;}
.veSpecifier .veTable .scroller .variableEditorRenderers {padding-right: 3px; -webkit-user-select: none; -moz-user-select: none; -ms-user-select: none;}
.veVariableValueSummary {display: inline-block;}
.veSpecifier .veTable .topHeaderWrapper,.veSpecifier .veTable .bottomRowHeaderWrapper {visibility: hidden; z-index: 0;}
.veMetaSummary {font-style: italic;}
.veSpecifier .veTable .scroller {overflow: hidden;}
.veSpecifier .veTable:hover .scroller {overflow: auto;}
.veSpecifier .veVariableName,.veSpecifier .veDimensionFont {font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px;}
.veSpecifier .veVariableName {padding-top: 3px;}
.veSpecifier .veDimensionFont {font-style: italic; color: #9A9A9A;}
.veSpecifier .scroller::-webkit-scrollbar-track {background-color: white;}
.veSpecifier .scroller::-webkit-scrollbar-corner {background-color: white;}
.veSpecifier .veTable .topRowHeaderWrapper {border: none; background-color: #F8F9FA;}
.veSpecifier .mw_type_ListBase.showCellBorders,.veSpecifier .veTable .topHeaderWrapper,.veSpecifier .veTable .bottomRowHeaderWrapper,.veSpecifier .veTable .verticalScrollSpacer,.veSpecifier .veTable .horizontalScrollSpacer {border: none;}
.veSpecifier .veTable .dataScrollerNode {border: 1px solid #BFBFBF;}
.veSpecifier .veTable .columnHeaderNode,.veSpecifier .veTable .rowHeaderNode,.veSpecifier .veTable .dataBody {font-family: Arial; font-size: 13px;}
.veSpecifier .veTable .columnHeaderNode,.veSpecifier .veTable .rowHeaderNode {color: #7F7F7F;}
.veSpecifier .veTable .dataBody {color: #000000;}
.veSpecifier .veTable .columnHeaderNode .cell .drag,.veSpecifier .veTable .columnHeaderNode .cell .header,.veSpecifier .veTable .topHeaderWrapper,.veSpecifier .veTable .bottomRowHeaderWrapper {background-color: #F8F9FA;}
.veSpecifier .veTable .columnHeaderNode .cell .dragBorder {border-right: 1px solid #F8F9FA;}
.veSpecifier .veTable .rowHeaderNode .cell {padding-top: 3px; text-align: center; border-bottom: 1px solid #F8F9FA;}
.veSpecifier .veTable .dataBody .cell .plainText {text-align: right;}
.veSpecifier .veTable .dataBody .row .cell {border-bottom: 1px solid #E9E9E9; border-right: 1px solid #E9E9E9;}
.embeddedOutputsVariableElement {font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px; white-space: pre-wrap; word-wrap: break-word; min-height: 18px; max-height: 250px; overflow: auto;}
.variableElement {}
.embeddedOutputsVariableElement.inlineElement {}
.inlineElement .variableElement {}
.variableValue { width: 100% !important; }
.embeddedOutputsMatrixElement {min-height: 18px; box-sizing: border-box; font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px;}
.matrixElement .variableValue {white-space: pre; display: inline-block; vertical-align: top; overflow: hidden;}
.embeddedOutputsMatrixElement.inlineElement {}
.embeddedOutputsMatrixElement.inlineElement .topHeaderWrapper {display: none;}
.embeddedOutputsMatrixElement.inlineElement .veTable .body {padding-top: 0 !important; max-height: 100px;}
.embeddedOutputsMatrixElement.inlineElement .veVariableName {padding-top: 0;}
.inlineElement .matrixElement {max-height: 300px;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer {white-space: nowrap; padding-top: 10px; margin-bottom: 3px;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer .horizontalEllipsis.hide,.embeddedOutputsMatrixElement .matrixElement .verticalEllipsis.hide {display: none;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer .horizontalEllipsis {margin-bottom: -3px;}
.matrixElement { max-height: none !important; }
.dijitTooltipDialog .dijitTooltipContainer .dijitTooltipContents .alertPlugin-alertMessage {min-width: 12px; max-width: 400px; max-height: 300px; overflow: auto;}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar {width: 11px; height: 11px;}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-track {background-color: rgba(0, 0, 0, 0);}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-corner {background-color: rgba(0, 0, 0, 0);}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-thumb {border-radius: 7px; background-color: rgba(0, 0, 0, 0.1); border: 2px solid rgba(0, 0, 0, 0); background-clip: padding-box;}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-thumb:hover {background-color: rgba(0, 0, 0, 0.2);}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-thumb {background-color: rgba(0, 0, 0, 0);}
.dijitTooltipDialog .alertPlugin-alertMessage:hover::-webkit-scrollbar-thumb {background-color: rgba(0, 0, 0, 0.1);}
.dijitTooltipDialog .alertPlugin-alertMessage:hover::-webkit-scrollbar-thumb:hover {background-color: rgba(0, 0, 0, 0.2);}
.alertPlugin-alertLine {position: absolute; display: initial; width: 40px; text-align: right; cursor: text;}
.alertPlugin-onTextLine {visibility: hidden;}
.alertPlugin-hasTooltip .alertPlugin-warningImg,.alertPlugin-hasTooltip .alertPlugin-errorImg {cursor: pointer;}
.alertPlugin-isStale {-webkit-filter: opacity(0.4) grayscale(80%); filter: opacity(0.4) grayscale(80%);}
.alertPlugin-alertLine .alertPlugin-errorElement {display: inline-block; margin-right: 4px;}
.alertPlugin-errorImg {position: relative; display: inline-block; background-color: rgb(204, 55, 41); width: 14px; height: 14px; margin-top: 1px; -webkit-border-radius: 7px;}
.alertPlugin-alertLine .alertPlugin-warningElement {display: inline-block; margin-right: 3px;}
.alertPlugin-warningImg {position: relative; display: inline-block; width: 16px; height: 15px; overflow: hidden;}
.diagnosticMessage-wrapper {font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {color: rgb(255,100,0); text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {color: rgb(230,0,0); text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart {white-space: pre-wrap;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {white-space: pre;}
.embeddedOutputsWarningElement{min-height: 18px; max-height: 250px; overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsErrorElement {min-height: 18px; max-height: 250px; overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.alertPlugin-warningImg {background-image: url()}
.alertPlugin-errorImg {background-image: url()}
</style></head><body><div class = "content"><div class = 'SectionBlock containment active'><h1 class = "S1"><span class = "S2">Modelling coupled miscible fluid flow and transport processes </span></h1><p class = "S3"><span class = "S2">Fluid flow in subsurface aquifers and reservoirs may be undertaken under variable fluid properties. For instance, fluid viscosity changes dynamically as a function of local fluid concentration, temperature, or pressure conditions. These changes induce a tight dependence between the flow and solute/heat transport processes. Many applications belong to these category:</span></p><ul class = "S4"><li class = "S5"><span class = "S0">In the context of miscible Enhanced Oil Recovery a less viscous fluid (i;e. cold water, hot water, carbon dioxide, surfactants, etc) is often injected into the host formation to assist in the oil recovery process which involves a number of injection and production wells geographically placed following a given pattern. The oil volume in less permeable regions in the reservoir is best recovered when the injected fluid has lower viscosity leading to early breakthrough and higher recovery in production wells. </span></li><li class = "S5"><span class = "S0">Aquifer remediation technologies require optimal placement of a number of active and passive wells to decontaminate an initially polluted site by driving the concentration distribution below a national (or international) target regulation threshold. Characterization of the chemical composition and viscosity of injected aqueous solution plays an important role in this design. </span></li><li class = "S5"><span class = "S0">Pressure control in a deep geothermal reservoir requires very often the design of injection and production doublet wells. Extracted hot fluid is cooled in an external facility to produce heat or energy and reinjected back into the host formation at a significantly lower temperature. Since the water (oe brine) viscosity is a strong function of temperature the injected fluid has significant higher viscosity and therefore lower mobility thus retarding the heat front dissipation and enhancing the half-life of the doublet. </span></li></ul><p class = "S3"><span class = "S2">The goal of this tutorial is to learn how to model these processes. We will consider a simple two-dimensional example with a highly heterogeneous distribution of the permeability and porosity fields. These fields are extracted from the first layer of the SPE10 comparative solution project realization. </span></p><p class = "S3"><span class = "S2">An injection and production wells are placed at the lower left and upper right corners of the computational domain. Next, three simulations will be performed: (i) a reference simulation with an equal viscosity ratio (ratio of injected fluid viscosity to that of the resident fluid), (ii) a second simulation with a gradually increasing viscosity ratio, and (iii)  a third simulation with a gradually decreasing viscosity ratio.  </span></p></div><p class = "S0"></p><div class = 'SectionBlock containment'><h2 class = "S6"><span class = "S2">Data preparation and simulations setup </span></h2><p class = "S3"><span class = "S2">Lest's start by constructing the computational domain:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% Domain size along X, Y &amp; Z directions</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">dom = RecDomain([365.76,670.56,0.6096]);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Dx = dom.Lx; Dy = dom.Ly; Dz = dom.Lz;</span></p></div></div><p class = "S11"><span class = "S2">Set the number of cells along X, Y, and Z directions: </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.Nx = 60;   Grid.Ny = 220;   Grid.Nz = 1; </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Nx = Grid.Nx;   Ny = Grid.Ny;    Nz = Grid.Nz; </span></p></div></div><p class = "S11"><span class = "S2">Calculate grid sizes along each space direction:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.hx = (Dx/Nx);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.hy = (Dy/Ny);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.hz = (Dz/Nz);</span></p></div></div><p class = "S11"><span class = "S2">Set total number of gridblocks:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.N = Grid.Nx*Grid.Ny*Grid.Nz;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">N = Grid.N;</span></p></div></div><p class = "S11"><span class = "S2">Set Grid data: </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">Grid.V = (Dx/Nx)*(Dy/Ny)*(Dz/Nz);                        </span><span class = "S9">% Cell volumes</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">Grid.K = 1.4e-11.*ones(3,Nx,Ny,Nz);                      </span><span class = "S9">% Unit-Darcy permeability</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">Grid.compr = 4.4e-10.*ones(Grid.Nx,Grid.Ny,Grid.Nz);     </span><span class = "S9">% compressibility</span></p></div></div><p class = "S11"><span class = "S2">Load porosity data from a saved mat file:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">load </span><span class = "S13">'..\data\spe_layer1_phi.mat' phi_layer</span><span class = "S10">;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">phi = phi_layer; clear </span><span class = "S13">phi_layer</span><span class = "S10">;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">phi = max (phi(:), 1e-3); </span><span class = "S9">% to avoid zero porosity cells</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.por = reshape(phi',Nx,Ny,Nz); </span></p></div></div><p class = "S11"><span class = "S2">Load permeability data from a saved mat file: </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">load </span><span class = "S13">'..\data\spe_layer1_perm.mat' perm_layer</span><span class = "S10">;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.K = perm_layer*9.869232667160130e-13;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.K = reshape(Grid.K',3,Nx,Ny,Nz); </span></p></div></div><p class = "S11"><span class = "S2">Draw porosity and log10-permeability maps: </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">G = cartGrid(0:Dx/Nx:Dx, 0:Dy/Ny:Dy);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">figure;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">subplot(1,2,1);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">PlotCellData(G, phi, </span><span class = "S13">'EdgeColor'</span><span class = "S12">,</span><span class = "S13">'none'</span><span class = "S10">);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">%colormap(jet(16));</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">axis </span><span class = "S13">tight equal</span><span class = "S12">; colorbar, title(</span><span class = "S13">'Porosity field'</span><span class = "S10">);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">subplot(1,2,2);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">PlotCellData(G, reshape(log10(Grid.K(1,:,:,:)),N,1), </span><span class = "S13">'EdgeColor'</span><span class = "S12">,</span><span class = "S13">'none'</span><span class = "S10">);</span></p></div><div class = 'inlineWrapper outputs'><p class = "S8 lineNode"><span class = "S12">axis </span><span class = "S13">tight equal</span><span class = "S12">; colorbar, title(</span><span class = "S13">'log10(K_x) permeability field'</span><span class = "S10">);</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsFigure" style="max-height: 800px; width: 889px;"><div class="figureElement"><img class="figureImage" draggable="false" src=""></div></div></div></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">clear </span><span class = "S13">phi</span><span class = "S10">;</span></p></div></div><p class = "S11"><span class = "S2">Next, add injection / production wells and their respective flow rates:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% Cell-centered injection/production flow rates</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">flow_rate  = 100/3600;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Qw         = zeros(N,1);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Qw([1 N])  = [+flow_rate -flow_rate];</span></p></div></div><p class = "S11"><span class = "S2">Initialize the only fluid used in this model:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% fluid properties units: density (Kg/m^3), initial viscosity (Pa.s)</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">water = Fluid(</span><span class = "S13">'water'</span><span class = "S10">,[1000, 1e-3]);</span></p></div></div><p class = "S11"><span class = "S2">Initialize fluid pressure in the computational domain to 100 bars:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">P</span><span class = "S10"> = 100e5.*ones(Nx,Ny,Nz);</span></p></div></div><p class = "S11"><span class = "S2">Define time stepping parameters including total simulation time, T, the number of time steps, nt, and the constant time step, dt in seconds</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">day = 3600*24;                     </span><span class = "S9">% seconds/day</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">T   = 5*day;                       </span><span class = "S9">% 5 days</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">nt  = 30;                          </span><span class = "S9">% number of time steps</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">dt  = T/nt;                        </span><span class = "S9">% time step</span></p></div></div><p class = "S11"><span class = "S2">Setup options for the Newton-Raphson solver:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">opt.tol     = 1e-9;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">opt.maxiter = 50;</span></p></div></div></div><p class = "S0"></p><div class = 'SectionBlock containment'><p class = "S3"><span class = "S2">We wil now solve a solute transport problem. A unit concentration is carried from the injection well at the lower left corner:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">D = 1e-9.*ones(3,Nx,Ny,Nz); </span><span class = "S9">% uniform diffusion coefficient through the domain </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">C1 = zeros(N,1); C1(1) = 1; </span><span class = "S9">% unit input concentration at the injector &amp; zero otherwise</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">Grid.sat = ones(Nx,Ny,Nz);  </span><span class = "S9">% unit saturation everywhere </span></p></div></div><p class = "S11"><span class = "S2">Go through the main diffusive-convective solute transport time loop. Advection and diffusion processes are simulated by splitting the advection and diffusive parts symmetrically as suggested by Strang to increase the concentration accuray by minimizing the time-splitting error:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">tic;                        </span><span class = "S9">% start timing the time marching loop</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% call TPFA flow solver for single-phase flow (we assume steady state flow here)</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">[~,V] = Pressure(Grid,water,Qw);   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S16">for </span><span class = "S10">t=1:nt</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    fprintf(</span><span class = "S13">'\nSolving solute transport problem. Time = %f days\n'</span><span class = "S10">, t*dt/day);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    </span><span class = "S9">% call the fully implicit first-order upwind convective transport for a half time step</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    C_bar = ImplicitConcentration(Grid,C1,V,Qw,dt/2,opt); </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">      </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    </span><span class = "S9">% solve implicitly the diffusion problem for this time step</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    [C_bar,~] = Diffusion(Grid,D,C_bar,dt);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    </span><span class = "S9">% do advection step for the remaining half time step</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    C1 = ImplicitConcentration(Grid,C_bar,V,Qw,dt/2,opt); </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    </span><span class = "S9">% note: C_bar is the intermediate concentration after the first advection step </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    </span><span class = "S9">%             and the diffusive step</span></p></div><div class = 'inlineWrapper outputs'><p class = "S8 lineNode"><span class = "S17">end</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 0.166667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 0.333333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 0.500000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 0.666667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 0.833333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.166667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.333333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.500000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.666667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.833333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.166667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.333333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.500000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.666667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.833333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.166667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.333333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.500000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.666667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.833333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.166667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.333333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.500000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.666667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.833333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 5.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div></div></div><div class = 'inlineWrapper outputs'><p class = "S8 lineNode"><span class = "S10">toc;</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Elapsed time is 8.663168 seconds.</div></div></div></div></div><p class = "S11"><span class = "S2">For this problem the numerical solution is obtained quickly after only few seconds. Letes examone the shape of the solute plume after 5 days of continuous injection. </span></p></div><p class = "S0"></p><div class = 'SectionBlock containment'><p class = "S3"><span class = "S2">Plot the final concentration plume with 11 filled contour lines corresponding to iso-concentrations from 0 to 1 with a 0.1 spacing:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">figure; </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% plot solute concentration for this first simulation</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),</span><span class = "S17">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),</span><span class = "S17">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">         reshape(C1,Nx,Ny)',</span><span class = "S13">'LevelList'</span><span class = "S10">,0:0.1:1);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">colormap(jet(16));</span></p></div><div class = 'inlineWrapper outputs'><p class = "S8 lineNode"><span class = "S12">axis </span><span class = "S13">tight equal</span><span class = "S12">; colorbar, title(</span><span class = "S13">'Concentration isocontours: Case 1'</span><span class = "S10">);</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsFigure" style="max-height: 800px; width: 889px;"><div class="figureElement"><img class="figureImage" draggable="false" src=""></div></div></div></div></div></div><p class = "S0"></p><div class = 'SectionBlock containment'><h2 class = "S6"><span class = "S2">Case 2: miscible flow and transport - increasing viscosity ratio</span></h2><p class = "S3"><span class = "S2">Now, we will repeat the last simultion but the flow and solute transport equations are coupled by a favorable viscosity ratio which is assumed to dependent linearly on the solute concentration following the Doles-Jones correlation:</span></p><p class = "S18"><span style="vertical-align:-7"><img src="" width="132" height="22" /></span></p><p class = "S3"><span class = "S2">where:</span></p><p class = "S3"><span style="vertical-align:-3"><img src="" width="36.5" height="18" /></span><span class = "S2"> is a correction term which depends on the nature of the solute being injected.</span></p><p class = "S3"><span style="vertical-align:-7"><img src="" width="34" height="22" /></span><span class = "S2"> is the viscosity at the reference concentration </span><span style="vertical-align:-7"><img src="" width="14.5" height="22" /></span><span class = "S2">.</span></p><p class = "S3"><span class = "S2">We define this viscosity dependence on the concentration with the help of an anonymous function. In the MATLAB scripting language an anonymous function is expressed as an analytic function handle (the @) followed by declaration of the dependent variables (herein c vector) and then the analytic function expression. In this simulation we set </span><span style="vertical-align:-3"><img src="" width="10" height="18" /></span><span class = "S2"> to +10 such that the viscosity magnitude will increase at a maximum by one order of magnitude (i.e. at c=1). </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">visc = @(c) 1e-3*(1+10.*c);</span></p></div></div><p class = "S11"><span class = "S2">Now, let's do a coupled miscible flow and transport simulation time stepping loop. The key difference with the previous simulation is that the flow solver is repeatedly called at the beginning of each time step with updated cell-centered viscosities in each grid cell. This task is simply performed by calling our previously defined anonymous viscosity function and setting this new array to target fluid property. Next, the solute transport algorithm is identical to that discussed in the previous case.</span></p><p class = "S3"><span class = "S2">Concentrations for this second case will be stored in C2 array to enable comparaison with results of first simulation stored in C1 array.</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">C2 = zeros(N,1); C2(1) = 1;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">tic;              </span><span class = "S9">% start timing the time marching loop</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S16">for </span><span class = "S10">t=1:nt</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    fprintf(</span><span class = "S13">'\nUpdating the flow.\n'</span><span class = "S10">);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    </span><span class = "S9">% viscosity recalculation </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    water.viscosity = visc(C2); </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    </span><span class = "S9">% call TPFA flow solver for single-phase flow </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    [P,V] = Pressure(Grid,water,Qw);   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    fprintf(</span><span class = "S13">'\nSolving solute transport problem. Time = %f days\n'</span><span class = "S10">, t*dt/day);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    </span><span class = "S9">% call the fully implicit first-order upwind convective transport for a half time step</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    C_bar = ImplicitConcentration(Grid,C2,V,Qw,dt/2,opt); </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    </span><span class = "S9">% solve implicitly the diffusion problem for this time step</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    [C_bar,~] = Diffusion(Grid,D,C_bar,dt);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    </span><span class = "S9">% do advection step for the remaining half time step</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    C2 = ImplicitConcentration(Grid,C_bar,V,Qw,dt/2,opt);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper outputs'><p class = "S8 lineNode"><span class = "S17">end</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsWarningElement" data-width="889" style="width: 889px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType"><div class="diagnosticMessage-messagePart">Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  7.017226e-17.</div><div class="diagnosticMessage-stackPart"></div></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 0.166667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsWarningElement" data-width="889" style="width: 889px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType"><div class="diagnosticMessage-messagePart">Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  1.724073e-16.</div><div class="diagnosticMessage-stackPart"></div></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 0.333333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 0.500000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 0.666667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 0.833333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.166667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.333333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.500000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.666667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.833333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.166667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsWarningElement" data-width="889" style="width: 889px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType"><div class="diagnosticMessage-messagePart">Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  5.323851e-17.</div><div class="diagnosticMessage-stackPart"></div></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.333333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsWarningElement" data-width="889" style="width: 889px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType"><div class="diagnosticMessage-messagePart">Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  7.239497e-17.</div><div class="diagnosticMessage-stackPart"></div></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.500000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.666667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.833333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.166667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.333333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.500000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.666667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.833333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsWarningElement" data-width="889" style="width: 889px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType"><div class="diagnosticMessage-messagePart">Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  1.453285e-16.</div><div class="diagnosticMessage-stackPart"></div></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.166667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.333333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.500000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.666667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsWarningElement" data-width="889" style="width: 889px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType"><div class="diagnosticMessage-messagePart">Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  3.587041e-18.</div><div class="diagnosticMessage-stackPart"></div></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.833333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 5.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div></div></div><div class = 'inlineWrapper outputs'><p class = "S8 lineNode"><span class = "S10">toc;</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Elapsed time is 8.691799 seconds.</div></div></div></div></div></div><p class = "S0"></p><div class = 'SectionBlock containment'><p class = "S3"><span class = "S2">Let's plot the new results:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">figure; </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% plot solute concentration </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),</span><span class = "S17">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),</span><span class = "S17">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">         reshape(C2,Nx,Ny)',</span><span class = "S13">'LevelList'</span><span class = "S10">,0:0.1:1);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">colormap(jet(16));</span></p></div><div class = 'inlineWrapper outputs'><p class = "S8 lineNode"><span class = "S12">axis </span><span class = "S13">tight equal</span><span class = "S12">; colorbar, title(</span><span class = "S13">'Concentration isocontours: Case 2'</span><span class = "S10">);</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsFigure" style="max-height: 800px; width: 889px;"><div class="figureElement"><img class="figureImage" draggable="false" src=""></div></div></div></div></div><p class = "S11"><span class = "S2">When injecting a more viscous fluid the solute plume advancement is less restricted so as shown in the last figure the concentration fronts are more compressed and a higher concentration gradient is observed at the advancing front unlike in the first case where it is more diffusive in high permeability regions. However, since the fluid is less mobile solute diffusion is more pronounced in less permeable zones than in the first case where the fluid residence time is shorter. Thus there is less fluid mixing between waters in more and less permeable zones in the first case than in the second one. </span></p></div><p class = "S0"></p><div class = 'SectionBlock containment'><h2 class = "S6"><span class = "S2">Case 3: miscible flow and transport - decreasing viscosity ratio</span></h2><p class = "S3"><span class = "S2">Now we will consider a similar problem except that the injected fluid is less viscous. By setting </span><span style="vertical-align:-3"><img src="" width="10" height="18" /></span><span class = "S2"> to -10  viscosity is variable over one order in magnitude in the opposite direction when considering the second case. </span></p><p class = "S3"><span class = "S2">Concentrations for this third simulation are stored in C3 array.</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">visc = @(c) 1e-3*(1-10.*c);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">C3 = zeros(N,1); C3(1) = 1;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">tic;              </span><span class = "S9">% start timing the time marching loop</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S16">for </span><span class = "S10">t=1:nt</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    fprintf(</span><span class = "S13">'\nUpdating the flow.\n'</span><span class = "S10">);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    </span><span class = "S9">% viscosity recalculation </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    water.viscosity = visc(C3); </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    </span><span class = "S9">% call TPFA flow solver for single-phase flow </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    [P,V] = Pressure(Grid,water,Qw);   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    fprintf(</span><span class = "S13">'\nSolving solute transport problem. Time = %f days\n'</span><span class = "S10">, t*dt/day);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    </span><span class = "S9">% call the fully implicit first-order upwind convective transport for a half time step</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    C_bar = ImplicitConcentration(Grid,C3,V,Qw,dt/2,opt); </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    </span><span class = "S9">% solve implicitly the diffusion problem for this time step</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    [C_bar,~] = Diffusion(Grid,D,C_bar,dt);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">    </span><span class = "S9">% do advection step for the remaining half time step</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">    C3 = ImplicitConcentration(Grid,C_bar,V,Qw,dt/2,opt);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper outputs'><p class = "S8 lineNode"><span class = "S17">end</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsWarningElement" data-width="889" style="width: 889px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType"><div class="diagnosticMessage-messagePart">Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  5.852943e-17.</div><div class="diagnosticMessage-stackPart"></div></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 0.166667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 0.333333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsWarningElement" data-width="889" style="width: 889px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType"><div class="diagnosticMessage-messagePart">Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  5.678964e-17.</div><div class="diagnosticMessage-stackPart"></div></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 0.500000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 0.666667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsWarningElement" data-width="889" style="width: 889px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType"><div class="diagnosticMessage-messagePart">Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  4.258956e-17.</div><div class="diagnosticMessage-stackPart"></div></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 0.833333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.166667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.333333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.500000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsWarningElement" data-width="889" style="width: 889px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType"><div class="diagnosticMessage-messagePart">Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  1.365170e-16.</div><div class="diagnosticMessage-stackPart"></div></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.666667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 1.833333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.166667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.333333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.500000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.666667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 2.833333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.166667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.333333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.500000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.666667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 3.833333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.166667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.333333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.500000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.666667 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 4.833333 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Updating the flow.</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving solute transport problem. Time = 5.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps<br>Newton-Raphson iteration converged in 1 steps</div></div></div></div><div class = 'inlineWrapper outputs'><p class = "S8 lineNode"><span class = "S10">toc;</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Elapsed time is 10.769796 seconds.</div></div></div></div></div></div><p class = "S0"></p><div class = 'SectionBlock containment'><p class = "S3"><span class = "S2">Let's plot the new results:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">figure; </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% plot solute concentration </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),</span><span class = "S17">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),</span><span class = "S17">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S12">         reshape(C3,Nx,Ny)',</span><span class = "S13">'LevelList'</span><span class = "S10">,0:0.1:1);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">colormap(jet(16));</span></p></div><div class = 'inlineWrapper outputs'><p class = "S8 lineNode"><span class = "S12">axis </span><span class = "S13">tight equal</span><span class = "S12">; colorbar, title(</span><span class = "S13">'Concentration isocontours: Case 3'</span><span class = "S10">);</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsFigure" style="max-height: 800px; width: 889px;"><div class="figureElement"><img class="figureImage" draggable="false" src=""></div></div></div></div></div><p class = "S11"><span class = "S2">The most important outcome of this third simulation is that the first concentration fronts of the solute plume when injecting a less viscous fluid are spreading much rapidly. As such, solute breakthrough in the production well is expected to occur much earlier than in the previous cases. </span></p></div></div>
<!-- 
##### SOURCE BEGIN #####
%% Modelling coupled miscible fluid flow and transport processes 
% Fluid flow in subsurface aquifers and reservoirs may be undertaken under variable 
% fluid properties. For instance, fluid viscosity changes dynamically as a function 
% of local fluid concentration, temperature, or pressure conditions. These changes 
% induce a tight dependence between the flow and solute/heat transport processes. 
% Many applications belong to these category:
% 
% * In the context of miscible Enhanced Oil Recovery a less viscous fluid (i;e. 
% cold water, hot water, carbon dioxide, surfactants, etc) is often injected into 
% the host formation to assist in the oil recovery process which involves a number 
% of injection and production wells geographically placed following a given pattern. 
% The oil volume in less permeable regions in the reservoir is best recovered 
% when the injected fluid has lower viscosity leading to early breakthrough and 
% higher recovery in production wells. 
% * Aquifer remediation technologies require optimal placement of a number of 
% active and passive wells to decontaminate an initially polluted site by driving 
% the concentration distribution below a national (or international) target regulation 
% threshold. Characterization of the chemical composition and viscosity of injected 
% aqueous solution plays an important role in this design. 
% * Pressure control in a deep geothermal reservoir requires very often the 
% design of injection and production doublet wells. Extracted hot fluid is cooled 
% in an external facility to produce heat or energy and reinjected back into the 
% host formation at a significantly lower temperature. Since the water (oe brine) 
% viscosity is a strong function of temperature the injected fluid has significant 
% higher viscosity and therefore lower mobility thus retarding the heat front 
% dissipation and enhancing the half-life of the doublet. 
% 
% The goal of this tutorial is to learn how to model these processes. We 
% will consider a simple two-dimensional example with a highly heterogeneous distribution 
% of the permeability and porosity fields. These fields are extracted from the 
% first layer of the SPE10 comparative solution project realization. 
% 
% An injection and production wells are placed at the lower left and upper 
% right corners of the computational domain. Next, three simulations will be performed: 
% (i) a reference simulation with an equal viscosity ratio (ratio of injected 
% fluid viscosity to that of the resident fluid), (ii) a second simulation with 
% a gradually increasing viscosity ratio, and (iii)  a third simulation with a 
% gradually decreasing viscosity ratio.  
%% Data preparation and simulations setup 
% Lest's start by constructing the computational domain:

% Domain size along X, Y & Z directions
dom = RecDomain([365.76,670.56,0.6096]);
Dx = dom.Lx; Dy = dom.Ly; Dz = dom.Lz;
%% 
% Set the number of cells along X, Y, and Z directions: 

Grid.Nx = 60;   Grid.Ny = 220;   Grid.Nz = 1; 
Nx = Grid.Nx;   Ny = Grid.Ny;    Nz = Grid.Nz; 
%% 
% Calculate grid sizes along each space direction:

Grid.hx = (Dx/Nx);
Grid.hy = (Dy/Ny);
Grid.hz = (Dz/Nz);
%% 
% Set total number of gridblocks:

Grid.N = Grid.Nx*Grid.Ny*Grid.Nz;
N = Grid.N;
%% 
% Set Grid data: 

Grid.V = (Dx/Nx)*(Dy/Ny)*(Dz/Nz);                        % Cell volumes
Grid.K = 1.4e-11.*ones(3,Nx,Ny,Nz);                      % Unit-Darcy permeability
Grid.compr = 4.4e-10.*ones(Grid.Nx,Grid.Ny,Grid.Nz);     % compressibility
%% 
% Load porosity data from a saved mat file:

load '..\data\spe_layer1_phi.mat' phi_layer;
phi = phi_layer; clear phi_layer;
phi = max (phi(:), 1e-3); % to avoid zero porosity cells
Grid.por = reshape(phi',Nx,Ny,Nz); 
%% 
% Load permeability data from a saved mat file: 

load '..\data\spe_layer1_perm.mat' perm_layer;
Grid.K = perm_layer*9.869232667160130e-13;
Grid.K = reshape(Grid.K',3,Nx,Ny,Nz); 
%% 
% Draw porosity and log10-permeability maps: 

G = cartGrid(0:Dx/Nx:Dx, 0:Dy/Ny:Dy);

figure;

subplot(1,2,1);
PlotCellData(G, phi, 'EdgeColor','none');
%colormap(jet(16));
axis tight equal; colorbar, title('Porosity field');

subplot(1,2,2);
PlotCellData(G, reshape(log10(Grid.K(1,:,:,:)),N,1), 'EdgeColor','none');
axis tight equal; colorbar, title('log10(K_x) permeability field');
clear phi;
%% 
% Next, add injection / production wells and their respective flow rates:

% Cell-centered injection/production flow rates
flow_rate  = 100/3600;
Qw         = zeros(N,1);
Qw([1 N])  = [+flow_rate -flow_rate];
%% 
% Initialize the only fluid used in this model:

% fluid properties units: density (Kg/m^3), initial viscosity (Pa.s)
water = Fluid('water',[1000, 1e-3]);
%% 
% Initialize fluid pressure in the computational domain to 100 bars:

P = 100e5.*ones(Nx,Ny,Nz);
%% 
% Define time stepping parameters including total simulation time, T, the 
% number of time steps, nt, and the constant time step, dt in seconds

day = 3600*24;                     % seconds/day
T   = 5*day;                       % 5 days
nt  = 30;                          % number of time steps
dt  = T/nt;                        % time step
%% 
% Setup options for the Newton-Raphson solver:

opt.tol     = 1e-9;
opt.maxiter = 50;
%% 
% We wil now solve a solute transport problem. A unit concentration is carried 
% from the injection well at the lower left corner:

D = 1e-9.*ones(3,Nx,Ny,Nz); % uniform diffusion coefficient through the domain 
C1 = zeros(N,1); C1(1) = 1; % unit input concentration at the injector & zero otherwise
Grid.sat = ones(Nx,Ny,Nz);  % unit saturation everywhere 
%% 
% Go through the main diffusive-convective solute transport time loop. Advection 
% and diffusion processes are simulated by splitting the advection and diffusive 
% parts symmetrically as suggested by Strang to increase the concentration accuray 
% by minimizing the time-splitting error:

tic;                        % start timing the time marching loop
    
% call TPFA flow solver for single-phase flow (we assume steady state flow here)
[~,V] = Pressure(Grid,water,Qw);   

for t=1:nt
    
    fprintf('\nSolving solute transport problem. Time = %f days\n', t*dt/day);
   
    % call the fully implicit first-order upwind convective transport for a half time step
    C_bar = ImplicitConcentration(Grid,C1,V,Qw,dt/2,opt); 
      
    % solve implicitly the diffusion problem for this time step
    [C_bar,~] = Diffusion(Grid,D,C_bar,dt);
   
    % do advection step for the remaining half time step
    C1 = ImplicitConcentration(Grid,C_bar,V,Qw,dt/2,opt); 

    % note: C_bar is the intermediate concentration after the first advection step 
    %             and the diffusive step
end
toc;
%% 
% For this problem the numerical solution is obtained quickly after only 
% few seconds. Letes examone the shape of the solute plume after 5 days of continuous 
% injection. 
% 
% Plot the final concentration plume with 11 filled contour lines corresponding 
% to iso-concentrations from 0 to 1 with a 0.1 spacing:

figure; 

% plot solute concentration for this first simulation
contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),...
         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),...
         reshape(C1,Nx,Ny)','LevelList',0:0.1:1);
colormap(jet(16));
axis tight equal; colorbar, title('Concentration isocontours: Case 1');
%% Case 2: miscible flow and transport - increasing viscosity ratio
% Now, we will repeat the last simultion but the flow and solute transport equations 
% are coupled by a favorable viscosity ratio which is assumed to dependent linearly 
% on the solute concentration following the Doles-Jones correlation:
% 
% $$\mu(C) = \mu(C_0)(1+\lambda C)$$
% 
% where:
% 
% $\lambda > 0$ is a correction term which depends on the nature of the solute 
% being injected.
% 
% $\mu(C_0)$ is the viscosity at the reference concentration $C_0$.
% 
% We define this viscosity dependence on the concentration with the help 
% of an anonymous function. In the MATLAB scripting language an anonymous function 
% is expressed as an analytic function handle (the @) followed by declaration 
% of the dependent variables (herein c vector) and then the analytic function 
% expression. In this simulation we set $\lambda$ to +10 such that the viscosity 
% magnitude will increase at a maximum by one order of magnitude (i.e. at c=1). 

visc = @(c) 1e-3*(1+10.*c);
%% 
% Now, let's do a coupled miscible flow and transport simulation time stepping 
% loop. The key difference with the previous simulation is that the flow solver 
% is repeatedly called at the beginning of each time step with updated cell-centered 
% viscosities in each grid cell. This task is simply performed by calling our 
% previously defined anonymous viscosity function and setting this new array to 
% target fluid property. Next, the solute transport algorithm is identical to 
% that discussed in the previous case.
% 
% Concentrations for this second case will be stored in C2 array to enable 
% comparaison with results of first simulation stored in C1 array.

C2 = zeros(N,1); C2(1) = 1;

tic;              % start timing the time marching loop
for t=1:nt
    
    fprintf('\nUpdating the flow.\n');
    
    % viscosity recalculation 
    water.viscosity = visc(C2); 
   
    % call TPFA flow solver for single-phase flow 
    [P,V] = Pressure(Grid,water,Qw);   
    
    fprintf('\nSolving solute transport problem. Time = %f days\n', t*dt/day);
   
    % call the fully implicit first-order upwind convective transport for a half time step
    C_bar = ImplicitConcentration(Grid,C2,V,Qw,dt/2,opt); 
   
    % solve implicitly the diffusion problem for this time step
    [C_bar,~] = Diffusion(Grid,D,C_bar,dt);
   
    % do advection step for the remaining half time step
    C2 = ImplicitConcentration(Grid,C_bar,V,Qw,dt/2,opt);

end
toc;
%% 
% Let's plot the new results:

figure; 

% plot solute concentration 
contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),...
         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),...
         reshape(C2,Nx,Ny)','LevelList',0:0.1:1);
colormap(jet(16));
axis tight equal; colorbar, title('Concentration isocontours: Case 2');
%% 
% When injecting a more viscous fluid the solute plume advancement is less 
% restricted so as shown in the last figure the concentration fronts are more 
% compressed and a higher concentration gradient is observed at the advancing 
% front unlike in the first case where it is more diffusive in high permeability 
% regions. However, since the fluid is less mobile solute diffusion is more pronounced 
% in less permeable zones than in the first case where the fluid residence time 
% is shorter. Thus there is less fluid mixing between waters in more and less 
% permeable zones in the first case than in the second one. 
%% Case 3: miscible flow and transport - decreasing viscosity ratio
% Now we will consider a similar problem except that the injected fluid is less 
% viscous. By setting $\lambda$ to -10  viscosity is variable over one order in 
% magnitude in the opposite direction when considering the second case. 
% 
% Concentrations for this third simulation are stored in C3 array.

visc = @(c) 1e-3*(1-10.*c);
C3 = zeros(N,1); C3(1) = 1;

tic;              % start timing the time marching loop
for t=1:nt
    
    fprintf('\nUpdating the flow.\n');
    
    % viscosity recalculation 
    water.viscosity = visc(C3); 
   
    % call TPFA flow solver for single-phase flow 
    [P,V] = Pressure(Grid,water,Qw);   
    
    fprintf('\nSolving solute transport problem. Time = %f days\n', t*dt/day);
   
    % call the fully implicit first-order upwind convective transport for a half time step
    C_bar = ImplicitConcentration(Grid,C3,V,Qw,dt/2,opt); 
   
    % solve implicitly the diffusion problem for this time step
    [C_bar,~] = Diffusion(Grid,D,C_bar,dt);
   
    % do advection step for the remaining half time step
    C3 = ImplicitConcentration(Grid,C_bar,V,Qw,dt/2,opt);

end
toc;
%% 
% Let's plot the new results:

figure; 

% plot solute concentration 
contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),...
         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),...
         reshape(C3,Nx,Ny)','LevelList',0:0.1:1);
colormap(jet(16));
axis tight equal; colorbar, title('Concentration isocontours: Case 3');
%% 
% The most important outcome of this third simulation is that the first 
% concentration fronts of the solute plume when injecting a less viscous fluid 
% are spreading much rapidly. As such, solute breakthrough in the production well 
% is expected to occur much earlier than in the previous cases.
##### SOURCE END #####
--></body></html>